Frontiers in Systems Biology

Logical models are widely used to study regulatory and signaling systems, yet their reuse, annotation, and exchange across tools remain challenging. Although SBML Level 3 Qualitative Models (SBML-qual) provides a standard representation, its XML-based syntax is difficult to inspect and edit directly. Here we introduce TabularQual, a spreadsheet-based, community-driven representation for Boolean and multi-valued logical models, together with a bidirectional converter between spreadsheets and SBML-qual. The converter is accessible programmatically and via a web interface to support diverse user workflows. We further describe two integration strategies that enable existing modeling tools to operate with TabularQual, either through SBML-qual exchange or via direct support. Case studies using the Stress Knowledge Map and CaSQ demonstrate how this integration supports model construction, curation, and reuse. TabularQual provides a practical bridge between human-readable model representations and standardized executable formats, supporting reproducibility, interoperability, and community-driven model development.

ObjectiveHigh-quality metadata is essential for scientific discovery, yet sparse annotations in rapidly growing repositories leave many biologically relevant details uncaptured. We evaluated whether large language models (LLMs) can accurately infer ion channel and receptor subtype metadata from source code in a neuroscience repository. Materials and MethodsWe extracted 5,133 model files from ModelDB. A subset of 1,100 was manually annotated; 253 were held out for testing, and the remainder split into training (80%) and validation (20%) sets. LLM-based approaches (GPT-5.2 and GPT-mini) were evaluated under zero-shot and heuristic-augmented prompting. Performance was assessed at type and subtype levels using accuracy, precision, recall, and F1 score. A feature-engineered XGBoost model using text- and simulation-derived features served as a baseline. ResultsLLMs outperformed the XGBoost baseline. At the type level, GPT-mini with heuristic augmentation achieved the highest performance (accuracy 96.0%, F1 0.962). At the subtype level, both GPT-5.2+heuristics and GPT-mini+heuristics achieved identical accuracy (88.1%), with GPT-5.2+heuristics achieving the highest F1(0.878). Model outputs were consistent across runs and errors confined to related mechanistic families. Discussion and ConclusionLLMs demonstrate strong potential for metadata annotation directly from source code, outperforming feature-engineering approaches with minimal tuning. However, performance varied across subtypes, and errors often reflected ambiguity or bias toward more common labels. These findings suggest LLMs may serve as practical tools for scalable metadata generation in biomedical repositories, although careful evaluation and domain-specific validation remain important. While demonstrated in computational neuroscience, this approach may generalize to repository-agnostic metadata annotation in other scientific code repositories.

Mechanistic ordinary differential equation models are widely used in systems biology to represent biochemical networks, population dynamics, cell-state transitions, and other biological processes; however, their predictive value depends critically on accurate parameter estimation from noisy and often sparse experimental data. In this tutorial, we present the Weak-form Estimation of Nonlinear Dynamics (WENDy) method as a forward-solver-free approach that reformulates parameter estimation as a covariance-corrected weak-form regression problem by integrating the model equations against compactly supported test functions. We present the background on the methodology through the lens of the familiar logistic equation, and we demonstrate applications of the method on real experimental data through two systems biology examples: a glycolytic oscillator with relatively dense time-course data and a sparse epithelial-mesenchymal cellstate transition model with multiple experimental replicates. Ultimately, using WENDy, we estimate interpretable biological parameters with uncertainty for systems with noisy and sometimes sparse available experimental data.

Endometrial receptivity assessment currently requires invasive tissue biopsy, yet recent randomized trials have questioned the clinical utility of biopsy-based approaches. Here we present EndoTwin-W, a four-layer mechanistic computational model that simulates human endometrial remodeling from hormone inputs through receptor binding, pathway scoring, and continuous-time Markov chain cell-state transitions across 17 cell states. Transition rates were optimized against scRNA-seq and microarray data, then validated by 5-fold cross-validation on an independent bulk RNA-seq cohort (n=236 biopsies), achieving significant correlations for 16 of 17 cell states (mean Spearman r = 0.505) with benchmark dominance over three null models for 13 of 17 states. The model identifies glycodelin-A (PAEP) and CA-125 (MUC16) as mechanistically grounded candidate circulating biomarkers capturing two principal receptivity failure modes: inadequate decidualization and excessive inflammation. Hill-function prediction of serum glycodelin-A shows strong rank-order calibration (Spearman rho = 0.833, p = 0.010). Cross-condition held-out validation against 9 independent datasets (244 samples) achieves significant concordance in 5 of 9 datasets (median rho = 0.435). A cross-dataset receptivity index analysis across 18 GEO datasets (21 comparisons) demonstrates mean AUC = 0.599 with correct direction in 76% of analyses, including significant RNA-seq validation (AUC = 0.770, p = 0.003). The divergence between predicted and measured biomarker values defines a Progesterone Resistance Score quantifying decidualization deficit and inflammation burden. EndoTwin-W provides a mechanistic framework and candidate blood-based biomarkers for receptivity assessment; prospective paired serum-tissue validation is required before clinical use. Author SummaryAssessing whether the uterine lining is ready for embryo implantation usually requires an invasive biopsy that is costly and cannot be repeated every cycle. We built a computer model called EndoTwin-W that simulates how ovarian hormones reshape the endometrium through hormone receptors, intracellular signaling, and changing cell states across the menstrual cycle. When we tested the model against published gene-expression datasets from hundreds of patient samples, it matched known endometrial cell states in 16 of 17 categories. Our main finding is that two blood proteins, glycodelin-A and CA-125, may serve as non-invasive markers of receptivity. Glycodelin-A reflects decidualization; CA-125 reflects inflammation. When the models predictions disagree with measured blood levels, the mismatch defines a two-dimensional progesterone resistance score that may help explain why some patients do not respond to progesterone despite normal hormone levels. We provide an open research website (https://endotwin-w.com: mirror: https://endotwinw.com) for exploration, but prospective clinical studies are still needed before this approach could guide patient care.

Glucosinolates are an important group of specialized metabolites in the Brassicaceae family, playing a role as defensive compounds against biotic attackers. In response to biotic stress, plants upregulate glucosinolate biosynthesis in part by increasing the abundance of enzymes in the glucosinolate biosynthetic pathway. As an increase in enzyme abundance is often preceded by an increase in the corresponding mRNA levels, the dynamic changes in mRNA levels should capture the information required to infer how metabolite levels change over time. In order to test this hypothesis, a time series of experimental glucosinolate content data collected from Arabidopsis thaliana, exposed to either a mock or methyl jasmonate (MeJA) treatment, as a proxy for biotic stress, was combined with existing mRNA abundance data over time at the same developmental stage and treatment. We propose the GEEM model, a multilevel mechanistic ordinary differential equation (ODE) model, which goes from Gene expression to an enzyme level model, followed by a Michaelis Menten kinetics metabolite model, to simulate the dynamics of a segment of the indolic glucosinolate pathway. In order to constrain the GEEM model, three models were fit to experimental de novo specialized metabolite data, using different degrees of freedom by utilizing both a Gradient Boosted Tree model with a tested architecture to predict the kinetic constants, and augmenting these predictions with a literature review of the known Michaelis Menten kinetic constants from the glucosinolate pathway. Using Sequential Monte Carlo - Approximate Bayesian Computing to fit the GEEM model to the experimental data, we showed that given the mRNA levels and initial concentrations of metabolites, the changes in specialized metabolites over time and treatment can be modeled. Author SummaryWe study how plants adjust their natural chemical defenses over time when they are under attack from living organisms. In the mustard family, including the subject of our experiment Arabidopsis, one important group of defense chemicals is called glucosinolates. When Arabidopsis is under attack, certain gene pathways can be activated or deactivated, allowing the plant to modulate the amount of enzymes they produce, which in turn modulates the levels of these defensive chemicals. In this work, we combine measurements of gene activity and glucosinolate levels from Arabidopsis treated with a compound used in stress signal that mimics insect or pathogen attack. We then constructed a mathematical model that goes from gene activity, to amount of enzyme present, and ends with the amounts of specific glucosinolates over time. By fitting this model to experimental data, we show that it is possible to predict how glucosinolate levels change over time from the gene activity and initial glucosinolate levels. Our approach offers a way to connect gene expression datasets to real changes in plant defense chemistry, with potential applications in plant breeding and insight into how these pathways change due to stress.

Heterogeneous biomedical knowledge networks (hetnets) integrate disparate data types, drugs, genes, diseases, and pathways, across independent sources; Hetionet (https://het.io) is a widely used example. A standard approach for assessing connectivity significance is XSwap, which permutes the hetnet P times and fits a gamma-hurdle null model to the degree-weighted path count (DWPC), pooling permuted values across pairs with matching source and target degrees to increase the effective sample size. This permutation approach has been highly successful in practice, but it faces four practical constraints in large graphs: (1) a finite resolution for the smallest reportable p-values, (2) computational cost that grows prohibitive at path lengths L [≥] 4 or 5, (3) a variance model (Var {propto} {micro}2) that departs from the configuration-model form (1 +{kappa} ){micro}, and (4) O(P 10m L) runtime. To complement this approach, we present HetNetEX (Heterogeneous Network EXact inference), which computes the null DWPC distribution analytically from degree sequences using the configuration model in O(Ln) time. In simulations at P = 200 across L = 1-4, HetNetEX achieves Spearman{rho} > 0.96 concordance with XSwap rankings while being >10,000x faster and providing analytical p-values without a resolution ceiling. High-degree pairs show larger XSwap sampling error than low-degree pairs, reflecting the finite-sample nature of permutation that analytical computation avoids.

The discovery of conserved regulatory motifs across different species is a fundamental challenge in systems biology, especially considering the noisy and incomplete nature of predicted gene regulatory networks (GRNs) and the intractability of the underlying graph alignment problem. Traditional network alignment methods frequently enforce one-to-one node mappings or strict topological isomorphism, which fail to accommodate the many-to-many orthology mappings caused by evolutionary gene duplication. Consequently, strict constraints often yield highly fragmented topological islands rather than cohesive functional modules. In this work, we propose a relaxed topological alignment algorithm designed to extract conserved regulatory structures from cross-species GRNs. We formulate the discovery process as a multi-objective optimization problem that balances sequence homology, functional coherence, and a normalized topological consensus. To navigate the exponentially scaling search space, we introduce a greedy seed-and-extend heuristic bounded by a dynamic{epsilon} -stopping condition, which evaluates marginal objective gains to prevent functional dilution. We validate our algorithm using time-series transcriptomic data from Arabidopsis thaliana, Zea mays, and Sorghum bicolor focused on drought and developmental stress responses. While a strict topological baseline extracted only fragmented subgraphs limited to 51 homologous tuples, our relaxed heuristic successfully converged on a highly connected 444-tuple module. The resulting topology effectively links strictly conserved upstream transcription factors to their highly duplicated, species-specific downstream pathways. Our algorithm provides a robust, scalable computational methodology for identifying core regulatory logic across complex biological systems, facilitating the translation of conserved network architectures among multiple species. Author summaryIdentifying shared regulatory mechanisms across diverse species is essential for understanding how complex biological systems evolve and adapt. However, traditional computer algorithms struggle to align these biological networks because evolution frequently duplicates genes, breaking simple one-to-one comparisons and producing highly fragmented results. To overcome this limitation, we developed a relaxed cross-species network alignment algorithm. Instead of demanding perfectly identical network shapes, our approach dynamically balances genetic sequence similarity, network structure, and biological function. We demonstrated the performance of our algorithm using plant drought-stress networks as a case study. While strict methods only found tiny, disconnected network fragments, our algorithm uncovered a functionally coherent, interconnected regulatory module across three distinct species. We discovered that while upstream command genes remain strictly conserved, they regulate highly customized, species-specific execution pathways downstream. Ultimately, our framework provides a scalable, species-agnostic method to decode complex systems, allowing researchers to translate conserved biological logic across diverse genomes.

The detection and monitoring of bat species using non-invasive sampling and molecular techniques has become increasingly popular in recent years. In Europe, these approaches have been applied to identify horseshoe bats of the genus Rhinolophus, which includes five species: R. hipposideros, R. ferrumequinum, R. euryale, R. mehelyi and R. blasii. While species-specific real-time PCR assays exist for R. ferrumequinum and R. hipposideros, no unified panel of real-time PCR assays currently enables the identification of all five European Rhinolophus species from non-invasively collected samples. Here, we developed five species-specific real-time PCR assays, each targeting interspecies nucleotide variation within the mitochondrial cytochrome b gene. To enhance single-base discrimination, RNase H-dependent PCR (rhPCR) primers were employed, incorporating cleavable blocked primers that require perfect complementarity for extension. The assays were applied to droppings non-invasively collected from 18 caves and one church in Serbia and Romania. Of the 149 samples analysed, 131 (88%) yielded successful amplification of Rhinolophus DNA. Detection probabilities for the three species identified in the field ranged from 0.49 to 0.82. Occupancy estimates varied, with R. euryale showing the highest (0.86; UI: 0.69-0.97) and R. mehelyi the lowest (0.23; UI: 0.08-0.43). The assays were capable of detecting up to three species concurrently within a single pooled sample (approximately 15 droppings). These assays are especially valuable for detecting R. mehelyi, given its rarity and uncertain distribution, and offer a robust tool for monitoring Rhinolophus populations across Europe.

Approximate Bayesian computation sequential Monte Carlo (ABC-SMC) propagates its particles with a perturbation kernel, and with the standard Normal kernel it degrades sharply as the parameter dimension grows, a failure usually attributed to dimension itself. We show instead that it is governed by the quality of the summary statistics, with dimension entering only through a separate and milder mechanism, and that the two must act together for the Normal kernel to break. The first ingredient is covariance overinflation: the kernel covariance, estimated from the particle cloud, overshoots the true posterior covariance by a factor set by information loss in the summary statistics. We derive this overscaling factor in closed form for a Gaussian model with sufficient statistics and show that it stays modest at any dimension, shrinking toward its baseline value as the tolerance tightens; the extreme values seen in practice (of order 103) are a signature of insufficient summaries, not of dimension. The second ingredient is perturbation overconcentration: the normalised Normal step size concentrates around one as the dimension grows, so every proposal overshoots by the same factor. Either ingredient alone is harmless; only their combination breaks the Normal kernel. A Cauchy kernel (multivariate t with one degree of freedom) removes the concentration, keeping a positive acceptance rate under arbitrary overscaling at a bounded worst-case cost of 1.87x in expected squared jump distance. In a Metropolis-Hastings framework we derive closed-form acceptance rates for both kernels that illustrate the advantage of the Cauchy kernel in this limit. A series of full ABC-SMC computational experiments on five problems at d = 12, including a hierarchical gene-expression model, show the Cauchy reducing the sliced Wasserstein distance to the reference posterior by factors of up to 50 with the same simulation budget. Since the summary statistics are commonly insufficient for the models that require ABC, overinflation is structural and the Cauchy perturbation kernel is the right default for problems in higher dimensions.

The aryl hydrocarbon receptor (AhR) is a ligand-activated transcription factor involved in xenobiotic sensing, as well as development, immunity, and tissue homeostasis. AhR signaling can proceed through a canonical and non-canonical pathway; the present study focuses on the canonical pathway. While ligand-dependent differences in binding affinities and direct ligand degradation kinetics are well known, and subtle differences in ligand binding can shape downstream signaling, it is still unclear which biochemical reaction steps within the canonical pathway are responsible for distinct ligand-specific transcriptional responses. Here, we developed a mechanistic ordinary differential equation model of the canonical AhR pathway. We calibrated the model to time-resolved qPCR measurements of Cyp1a1 and Ahrr mRNA in mouse bone-marrow-derived macrophages exposed to structurally diverse, environmentally relevant ligands with known immunomodulatory activity (3-methylcholanthrene, indolo[3,2-b]carbazole, and bisphenol A) using global optimization under a heteroskedastic likelihood. To dissect ligand specificity, we evaluated 528 candidate models that allow one or two ligand-involving reaction rate constants to vary. Akaike-based model selection reveals a dominant dynamical regime governed by promoter occupancy and target-gene mRNA synthesis, indicating that ligand-specific transcriptional responses are primarily encoded at the level of transcriptional regulation rather than upstream signaling events. The resulting model is made available in SBML and PEtab formats for reproducibility, and to enable further research into whether ligand-specific effects are conserved or rewired across cell types.

AO_SCPLOWBSTRACTC_SCPLOWThe stochastic simulation algorithm (SSA) is widely used to perform exact forward simulation of discrete stochastic processes in biology. However, the computational cost, driven by sequential event-by-event sampling across large ensembles, remains a computational barrier. We investigate whether reduced-precision floating-point arithmetic can accelerate SSA without degrading statistical fidelity, drawing on the success of reduced-precision methods in weather and climate modelling. We evaluate two strategies across five canonical models (birth-death, Schlogl, Telegraph, dimerisation, repressilator): (i) mixed precision, computing propensities in 16-bit while maintaining accumulators in 32-bit; and (ii) uniform precision, performing all arithmetic in 16-bit. Mixed-precision SSA produces ensemble statistics that closely match the 64-bit reference for all models, as measured by Kolmogorov-Smirnov tests and Wasserstein distances. Under uniform precision, deterministic rounding introduces systematic biases across several models, with catastrophic failures in some cases. Stochastic rounding (SR) and propensity normalisation eliminate these biases, restoring distributional fidelity across all models tested (KS p > 0.05). Our results establish mixed-precision SSA with SR as a viable acceleration strategy for mathematical biology: 16-bit formats shrink per-variable data size by 2-4x relative to fp32/fp64, yielding comparable reductions in memory footprint and up to ~ 1.5x wall-clock speedup on CPU hardware that lacks native 16-bit arithmetic. As a hardware-level acceleration, mixed-precision SSA complements algorithmic methods such as tau-leaping and maps naturally onto modern GPU and TPU architectures with native 16-bit arithmetic.

Understanding how gene regulatory networks respond to global cell perturbations remains a central challenge in systems biology and network inference. Modular Response Analysis (MRA) provides a mathematical framework to infer gene-to-gene directed connectivity graphs from perturbation experiments; however, classical MRA captures direct gene-to-gene influences, and does not explicitly account for global stimuli that simultaneously change the graph. Here, we introduce MRA+, an extension of MRA, that incorporates the effect of global perturbations into gene-to-gene graph inference. MRA+ assumes a sequential experimental design in which targeted gene perturbations are followed by the application of a global stimulus, enabling the separation of connectivity changes from direct gene induction. The method estimates network connectivity under induced conditions and quantifies gene-specific induction strengths, which represent contributions to expression changes arising from mechanisms external to the inferred network. In the case of single-cell expression data, we present a bootstrap strategy to assess the robustness of inferred connectivity coefficients and propose a complementary criterion based on sign stability to interpret weak or non-significant estimates. Together, these developments provide a general framework for robust inference of gene connectivity graphs in the presence of global perturbations, applicable to diverse biological and experimental contexts.

Cells communicate via extracellular ligands, such as hormones, which bind to plasma membrane receptors and trigger intracellular signaling cascades. G Protein-Coupled Receptors (GPCRs) exemplify this mechanism by initiating signaling both at the cell surface and, from intracellular compartments such as endosomes. The kinetics and spatial localization of these signals are critical determinants of cellular responses, yet receptor trafficking-including internalization, endosomal sorting, and recycling-remains a pivotal but often overlooked component of theoretical GPCR models. In this study, we present a mathematical framework that integrates receptor trafficking and signaling compartmentalization into generic GPCR dynamic models. Using a compartmentalized approach based on systems of ordinary differential equations (Chemical Reaction Networks), we analyze how receptor internalization and recycling modulate ligand-induced responses. Our results show that the balance between plasma membrane and endosomal signaling can significantly enhance or diminish ligand efficacy. Calibrated with high-throughput kinetic data, our model offers a refined tool for ligand pharmacological characterization and advances the understanding of GPCR signaling spatial organization.

BackgroundFederated learning (FL) enables collaborative model training across institutions without sharing patient-level data. However, standard FL algorithms such as FedAvg degrade under non-independently and non-identically distributed (non-IID) data, a prevalent condition when patient demographics, scanner hardware, and disease prevalence differ across hospital sites. ObjectiveWe propose iPS-MFFL (Individualized Per-Site Meta-Federated Feature Learning), a federated framework with a hierarchical local-model architecture that addresses non-IID heterogeneity through (1) a shared feature extractor, (2) multiple weak-learner classification heads that can be trained with heterogeneous training objectives to promote complementary decision boundaries, (3) independent per-learner server aggregation so that each weak learners parameters are averaged only with its counterparts at other clients, and (4) a lightweight meta-model -- itself federated -- that adaptively stacks the weak-learner outputs. The specific choices of backbone, weak-learner training objectives, and meta-model are implementation details; in this work we use an ImageNet-pretrained ResNet18 and three heterogeneous losses as a concrete instantiation. MethodsWe evaluate on the Brain Tumor MRI Classification dataset (7,200 images; 4 classes: glioma, meningioma, pituitary tumor, no tumor) partitioned across K = 5 simulated hospital sites using Dirichlet non-IID sampling ( = 0.3). Four baselines are compared: Local-only training, FedAvg, FedProx, and Freeze-FT. All experiments are repeated over three random seeds (13, 42, 2025) and evaluated using paired t-tests, Cohens d effect sizes, and post-hoc power analysis. ResultsiPS-MFFL achieved the highest mean final-round test accuracy point estimate of 85.42 {+/-} 8.70% (mean {+/-} SD across three seeds), compared to FedAvg (78.48 {+/-} 12.66%), FedProx (78.33 {+/-} 14.64%), Freeze-FT (73.98 {+/-} 21.09%), and Local (58.10 {+/-} 11.77%). iPS-MFFL showed the smallest cross-seed SD, suggesting greater robustness to partition heterogeneity. However, one-way ANOVA did not reach statistical significance (F = 1.52, p = 0.270), reflecting the limited number of experimental seeds. Cohens d effect sizes relative to iPS-MFFL ranged from 0.59 (vs. FedProx) to 2.64 (vs. Local); post-hoc pairwise comparisons are reported as exploratory given the non-significant omnibus test. Post-hoc power analysis indicated that statistical power for FL baseline comparisons was only 0.10-0.12 for the observed effect sizes (d {approx} 0.6) at n = 3 seeds. ConclusionsiPS-MFFL provides a practical approach to heterogeneous federated brain tumor classification by combining transfer learning, contrastive weak-learner diversity, and meta-learning. The framework demonstrated the highest mean accuracy and lowest variance across diverse data partitions. Validation with larger seed pools ([≥] 10 seeds for 80% power), ablation studies, and external multi-center cohorts is needed to establish generality.

Pairwise analysis of genomic and spatial data offers opportunities to identify and estimate the associations between covariates and the transmission of pathogens between individuals. However, such pairwise analyses are computationally intensive, and may not be feasible to conduct given the high dyad count in even moderately sized datasets. Here we compare two approaches to increase the efficiency of pairwise analysis for large datasets. We quantify and compare the performance of divide-and-conquer Bayesian model fitting and pairwise case-control approaches for estimating associations between individual- and pair-level covariates and shared membership in a transmission cluster. We utilize a large dataset (n=4,154) of spatially-referenced, genomically-sequenced Mycobacterium tuberculosis isolates collected from a single city for this analysis. We find that the case-control approach produces unbiased estimates of effect sizes with expected credible interval coverage and is more robust than the divide-and-conquer method when effect sizes are large. Thus, we recommend using the case-control approach with at least three controls per case to downscale datasets for pairwise analysis when analysis of the entire dataset is not possible. This approach mitigates the computational challenges of pairwise Bayesian modeling on datasets that require significant computational resources while maintaining desired inferential properties. Author SummaryPairwise analyses of large datasets to study pathogen transmission are computationally demanding because they typically require simultaneous analysis of each possible pair of individuals in a dataset; as datasets become larger these analyses often are not feasible to conduct even with access to high-performance computing resources. In this work, we compare a case-control approach and divide-and-conquer approaches for more efficient pairwise analysis of large datasets. Using a large dataset of Mycobacterium tuberculosis isolates including genetic and spatial data, we investigate the performance of each method for estimating the associations between host covariates and genetic clustering of isolates. We find that the case-control approach is generally preferred over methods which first divide the data into subsets and then combine results. While additional extensions of these analyses are needed to test the generality of these findings to other data settings, this work provides a practical way forward for the pairwise analysis of large datasets to study pathogen transmission.

The osteocyte network in bone is believed to play an important role for how bone tissues sense and respond to mechanical stimulation. Yet, bone adaptation to mechanical loads is often conceptualised as a simple response to mechanical stimuli, such as Wolffs law, which is based on mechanical variables only and takes no account of the cellular basis of mechanosensation. Wolffs law presumes the existence of a reference mechanical stimulus, the mechanical setpoint, above which bone is consolidated, and under which bone is removed. In this paper, we develop a theory of bone tissue sensing and adaptation based on osteocytes to provide new understanding of the role played by osteocyte signals in mechanical adaptation. In this theory, the mechanical setpoint of Frosts mechanostat is explicitly embodied as osteocyte properties involved in mechanotransduction. The mechanical setpoint is allowed to adapt due to the replacement of osteocytes during remodelling, making the setpoint space and time dependent. We propose a mathematical model to implement this new theory of bone adapation and present numerical simulations of this model to explore how mechanobiological response curves (effective Wolffs laws) are modulated by setpoint adaptation during remodelling. By accounting for varying osteocyte populations within bone tissue, we explore bone adaptation under osteocyte disruptions, which is particularly relevant to age-related bone loss. Our model suggests that biological disruptions of remodelling balance cannot always be compensated by mechanical feedback, and that setpoint adaptation during remodelling may have significant observable consequences, such as hysteresis in bone response signatures that resemble lazy zones.

1Endogenous circadian clocks control circadian rhythms in physiology and behavior. The predominant hypothesis of biological timing suggests that the responsible master clock for all endogenous circadian rhythms is constituted by an evolutionary conserved transcriptional-translational feedback loop (TTFL) clock consisting of positive feedforward and negative feedback elements. Unexpectedly, in contrast to the evolutionary derived insect Drosophila, RNAi-dependent knockdown of any of the negative feedback elements of the core TTFL clock in the basal Madeira cockroach Rhyparobia maderae does not delete circadian rhythms in locomotor activity. Shown here, neither RNAi-dependent triple knockdowns of all three negative feedback elements Period, Timeless 1, and Cryptochrome 2, nor single and double knockdown of the positive elements Clock and Cycle did directly delete circadian locomotor rhythms as mRNA levels declined. Thus, our experimental data do not support the predominant hierarchical hypothesis of circadian timing. To explore alternative mechanisms, we constructed a computational model of a neuronal circadian pacemaker network using planar switching affine systems (PSAS). The PSAS model comprises plasma membrane-associated posttranslational feedback loop (PTFL) clocks that are coupled to the TTFL nuclear clocks. Modeling results aligned with our experimental results. Therefore, both our experimental and modeling data support a systemic hypothesis of biological timing. 3 Significance statementBased mostly upon genetic studies in derived taxa like Drosophila it is hypothesized that circadian timing of behavior is strictly controlled by specific circadian clock neurons in the brain, realized through a transcriptional-translational feedback loop (TTFL) clock. In contrast to this common hierarchical model that requires transcription, we provide evidence in a basal taxon - the Madeira cockroach - for a systemic explanation of circadian timing of behavior that is based on coupled TTFL and posttranslational feedback loop (PTFL) clocks in adaptive neuronal networks.

We present the BCCWJ-Brain dataset, a multi-modal neuroimaging resource comprising functional magnetic resonance imaging (fMRI), magnetoencephalography (MEG), and electroencephalography (EEG) data recorded from native Japanese speakers reading newspaper articles from the Balanced Corpus of Contemporary Written Japanese (BCCWJ). Neural data were collected from 112 participants (36 fMRI, 35 MEG, and 41 EEG) as they read twenty newspaper articles presented in a Rapid Serial Visual Presentation (RSVP) paradigm. By providing three complementary neuroimaging modalities collected under identical naturalistic reading stimuli, this dataset provides a cognitive benchmark for computational models such as large language models. The dataset is publicly available on the OpenNeuro platform, offering a valuable resource for neuroscience, natural language processing, and related research fields.

Many biological characteristics arise by interactions between more than one biological organism or unit. Fertilization success in sexually reproducing species represents such an extended phenotype where both mates are required to be fertile for a successful outcome. Consequently, predictive models should account for the joint nature of reproductive performance while offering interpretable estimates for individual mate contributions. Recent advances in genomics and machine learning (ML) provide standardized, high-dimensional genetic information on one hand and computational tools capable of modeling complex biological systems on the other. Here, we construct and evaluate two-tower (TT) machine learning architectures for genomic prediction of binary reproductive outcomes and recovery of sex-specific fertility liabilities. Simulated datasets, generated under a range of genetic architectures, were utilized to compare multilayer perceptron (TT-MLP), convolutional neural network (TT-CNN), and L1-regularized linear (TT-LASSO) two-tower models. Simulation scenarios varied sex-specific heritabilities, genetic correlations, infertility prevalence, mating structure, and sex-specific infertility rates. Models were evaluated with regard to their ability to predict reproductive success at pair level and also recover true underlying genetic values for male and female fertility. Prediction accuracy increased with the underlying heritable component as expected, while sex-specific tower-scores successfully recovered latent fertility liabilities despite models being trained only on observed joint outcomes. TT-LASSO achieved the highest overall classification performance, whereas TT-MLP provided more balanced and consistent recovery of sex-specific genetic values across scenarios. An additional simulation, incorporating genotype-dependent mate compatibility demonstrated advantages of fully-connected neural networks for capturing non-additive interactions. These results indicate that two-tower frameworks provide a powerful approach for modeling reproductive traits, enabling simultaneous prediction of aggregate reproductive outcomes and sex-specific fertility liabilities from genotypic information.

Quantitative Systems Pharmacology (QSP) models play an inherently interventional role in pharmaceutical research and development, functioning as executable causal systems for designing, evaluating, and replacing clinical trials. However, deploying QSP as an experimental planning engine remains constrained by the difficulty of translating unstructured literature descriptions of clinical or preclinical scenarios into reproducible, simulation-ready model interventions. Motivated by this issue, we propose an agent-based framework that operationalizes QSP models as intervention-ready experimental systems by automatically extracting and executing literature-derived scenarios. The framework combines semantic grounding of model entities with a large language model (LLM)-driven Scenario Extractor and a dual-agent Scenario Mapper. Rather than relying on opaque, single-shot reasoning, our pipeline converts free-text interventions into precise parameter configurations through discrete, verifiable work orders. Moreover, our dynamic Human-in-the-Loop (HITL) strategy empowers modelers to resolve biological ambiguities interactively. Across four diverse kinetic ordinary differential equation (ODE)/QSP models and seven Subject Matter Expert (SME)-curated literature scenarios, our model resolved all selected scenarios into correct executable parameter changes, including multi-dose interventions, unit conversions, no-op scenarios, and ambiguity-triggered HITL cases, demonstrating that structured collaboration between experts and agentic systems can resolve scenarios that standalone raw Systems Biology Markup Language (SBML) reasoning LLM calls handle unreliably.